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Abstract 

We study the generation of gravitational waves in the electroweak phase tran- 
sition. We consider a few extensions of the Standard Model, namely, the addition 
of scalar singlets, the minimal super symmetric extension, and the addition of TeV 
fermions. For each model we consider the complete dynamics of the phase transition. 
In particular, we estimate the friction force acting on bubble walls, and we take into 
account the fact that they can propagate either as detonations or as deflagrations 
preceded by shock fronts, or they can run away. We compute the peak frequency 
and peak intensity of the gravitational radiation generated by bubble collisions and 
turbulence. We discuss the detectability by proposed spaceborne detectors. For 
the models we considered, runaway walls require significant fine tuning of the pa- 
rameters, and the gravitational wave signal from bubble collisions is generally much 
weaker than that from turbulence. Although the predicted signal is in most cases 
rather low for the sensitivity of LISA, models with strongly coupled extra scalars 
reach this sensitivity for frequencies / ~ 10~ 4 Hz, and give intensities as high as 

h 2 fl GW ~ io- 8 . 
1 Introduction 

Several gravitational wave (GW) detectors are currently planned to be constructed in 
space [1, 2, 3, 4, 5, 6]. The laser interferometer space antenna (LISA) [2, 3] is designed to 
detect the passage of a gravitational wave by measuring the time- varying changes of optical 
pathlength between free-falling masses. LISA consists of three spacecraft in heliocentric 
orbits, forming a triangle with sides ~ 10 9 m long. The LISA program was born more 
than ten years ago as a joint project of ESA and NASA. Recently, a variant of LISA 
was proposed, which is called New Gravitational wave Observatory (NGO) or evolved 
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LISA (eLISA). The Big Bang observer (BBO) [4] has been proposed as the successor 
of LISA. BBO is composed of four LISA type space detectors orbiting the sun, two of 
them collocated. In this case the arm length is ~ 10 7 m. A Japanese project with 
similar characteristics is the Deci-Hertz Interferometer Gravitational-wave Observatory 
(DECIGO) [5, 6]. The latter detectors would bridge the frequency gap between LISA and 
terrestrial detectors. 

These GW observatories will be able to measure a stochastic background of cosmolog- 
ical origin [7, 8, 9]. The detection of a primordial background of gravitational radiation 
would provide a direct probe of the physics in the early Universe, since GWs propa- 
gate freely after being produced. Cosmological sources of gravitational radiation include 
quantum fluctuations during inflation (see, e.g., [10]), scalar condensate fragmentation 
into Q-balls [11], cosmic strings (see, e.g., [12]), plasma turbulence and magnetic fields 
(see, e.g., [13, 14]). A possible scenario for the generation of a primordial GW background 
is a first-order phase transition of the Universe [15, 16, 17]. Quite interestingly, GWs pro- 
duced at the temperature scale of the electroweak phase transition, T* ~ 100 GeV , would 
have a characteristic frequency today (after redshifting) near the sensitivity peak of LISA, 
/ ~ lmHz. This motivated the investigation of GW production in the electroweak phase 
transition [18, 19, 20, 21, 22, 23]. 

In a first-order phase transition, bubbles of the stable phase nucleate and expand, 
converting the high-temperature phase into the low-temperature one (for the dynamics 
of a cosmological first-order phase transition see, e.g., [24, 25, 26] and references therein). 
Gravitational waves are generated either by the collisions of bubbles [15, 16, 17, 27, 28] or 
by the turbulence that is produced in the plasma due to the motion of bubble walls [17, 29, 
30, 31, 32]. In general, turbulence turns out to be a more effective source of gravitational 
radiation than bubble collisions. In the Standard Model (SM), the electroweak phase 
transition is not first-order [33]. As a consequence, the disturbance caused in the fluid 
is not enough to generate a significant GW signal. Models which give strongly first- 
order phase transitions and, consequently, a greater departure from equilibrium have 
been extensively studied in the context of electroweak baryogenesis [34]. 

The GW production can be calculated as a function of a few quantities related to 
the dynamics of the phase transition (see, e.g., [13, 14, 15, 16, 17, 27, 32]). These quan- 
tities are the temperature, the bubble wall velocity, the bubble size (or the duration of 
the phase transition), and the fraction of the released energy which goes into bulk mo- 
tions of the fluid (the efficiency factor). The latter can be calculated as a function of 
the amount of supercooling and the wall velocity [17]. Furthermore, the wall is usually 
assumed to propagate as a Jouguet detonation, so that the wall velocity and the efficiency 
factor have a simple dependence on the amount of supercooling. This motivated some 
model-independent analysis which find the electroweak GW spectrum as a function of 
two parameters, namely, the duration of the phase transition and the amount of super- 
cooling (see, e.g., [19]). In a given model, though, these parameters are linked, and it is 
important to investigate specific cases. Such investigations were performed, e.g., in Refs. 
[18, 20, 22, 23]. 

However, it is well known that the wall velocity does not only depend on the amount 
of supercooling but also on the friction with the surrounding plasma. In general, the 
hydrodynamic solution is not a Jouguet detonation. Contrary to the case of gravitational 
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waves, calculations of electroweak baryogenesis usually assume small wall velocities. The 
bubble growth mechanism has been studied for several years (see, e.g., [35]). Recently, 
there has been a renewed interest. The wall velocity was calculated taking into account 
hydrodynamics and microphysics in Refs. [36, 37, 38]. In Ref. [39] the microphysics 
in the ultra-relativistic regime was considered, finding that bubble walls may run away. 
The efficiency factor was calculated as a function of the wall velocity in Ref. [21] for 
deflagrations and, more recently, in Refs. [38, 40] for the whole range of wall velocities. 

In the present paper we study the generation of gravitational waves in the electroweak 
phase transition. We consider physical models and we include in the calculation some 
aspects of the dynamics which have not been taken into account previously. In particular, 
we incorporate the recent results on the hydrodynamics and microphysics of moving walls. 
We follow the evolution of the phase transition, taking into account the nucleation and 
expansion of bubbles, and the variation of temperature. We also take into account the 
effect of temperature inhomogeneities on the nucleation rate. We consider extensions of 
the SM with extra bosons, extra fermions, and the MSSM in the light-stop scenario. Our 
aim is to discuss the detectability of the gravitational radiation by LISA and other pro- 
posed detectors. Thus, we calculate the peak of the GW spectrum from bubble collisions 
and from turbulence, as a function of the parameters of each model. In this work we 
shall ignore the possible presence of magnetic fields, which would modify the turbulence 
mechanism [32, 41]. 

The plan of the paper is the following. In the next section we review the mechanisms 
for generation of gravitational waves by turbulence and bubble collisions. In section 3 we 
consider the nucleation, expansion and collisions of bubbles, and the energy injected into 
bulk motions of the fluid. We shall use results from Refs. [36, 37, 42] for the wall velocity, 
and results from Ref. [40] for the kinetic energy of the fluid. In section 4 we write down 
the one-loop finite-temperature effective potential which we shall use to calculate the 
thermodynamic parameters and the evolution of the phase transition. We also consider 
the general expression for the friction and the condition for "runaway" walls. In section 
5, we solve the equations for the dynamics of the phase transition and calculate the peak 
frequency and energy density of GWs. In section 6 we compare the signals obtained 
for the different models, and we discuss the possibility of observation at several planned 
space-based GW antennas. Finally, in section 7 we summarize our conclusions. 



The energy density of gravitational radiation is usually expressed in terms of the quantity 



where pew is the energy density of the GWs, / is the frequency, and p c is the critical 
energy density today, defined by p c = 3Hq/(8tcG), where H = 100 /ikms~ 1 Mpc~ 1 is 
the present day Hubble expansion rate, with h ~ 0.72 [43], and G is Newton's constant. 
Alternatively, the GW spectrum is often given in terms of the characteristic amplitude h c 
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or the root spectral density y/S. These are related to h 2 il GW by [1] 



^cw = *° ) , (2) 



2 

Hz 1.263 x 10- 18 , 
h c = y/2fS. (3) 



We shall consider the generation of GWs by bulk motions of the plasma during the 
electroweak phase transition. 

2.1 General features of GWs from bulk motions of the plasma 

For GWs originated at a time t*, a frequency /* redshifted to today is given by /o = 
/*a*/ao, where the ratio of the scale factor at t — t* to the scale factor today is given 
by the adiabatic expansion relation (goT^) / (g*T^) = al/a^, where go, To and g*,T* are 
the number of relativistic degrees of freedom (d.o.f.) and the temperature today and at 
t — respectively. We have 

a * o m-16 A00\ 1/3 100GeV 

— w 8 x 10 — — . (4) 

a V 9* J T * 

The typical wavelength will be a fraction of the Hubble size H^ 1 . Therefore, it is conve- 
nient to consider /*/#*, where the Hubble rate is given by the Friedmann equation, 

Hi = —p.. (5) 

Here, p* = p^ + p vac is the total energy density, where p^ = n 2 g*T^/ '30 is the energy 
density of radiation and p vac is the false vacuum energy density. Thus, we can express the 
frequency of the GWs today as 

The characteristic frequency is determined by the typical length scale of the source L s . 
Thus, one expects the peak of the spectrum to be at a frequency /* ~ l/As- 

The energy density of gravitational waves is given by [7] p G w( x ^) ~ (pth^dth^) / G , 
where is the tensor metric perturbation, and the brackets denote ensemble average. 
The equation for h^ u is of the form Oh^ u ~ GT^ U , where T^ u is the energy- momentum 
tensor of the source. On dimensional grounds, one expects the magnitude of h^ u to be 
given by Lg 2 h ~ Gpx, where px is the average kinetic energy density in bulk motions 
of the relativistic fluid. Similarly, we expect d t h ~ GpxLs- Therefore, we have pew ~ 
Gp\l? s . Using Eq. (5), this gives pew* ~ (pk/p*) 2 (LsH*) 2 p*. After the phase transition, 
the total energy density p* goes into radiation, and then evolves as (a*/ao) 4 . Since pew 
dilutes like radiation, we have today pcwo ~ (pk/p)1 (LsH)1 Pro- The energy density of 
radiation today is given by pm/pc = ~ 5 x 10~ 5 [43]. This estimate gives, for the 
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amplitude of the GW spectrum today 1 , 




In a first-order phase transition, the moving walls of expanding bubbles cause perturba- 
tions in the cosmic fluid. The bulk motions of the fluid produce GWs once bubbles collide 
and lose their spherical symmetry. In addition, bubble collisions generate turbulence, 
which is another source of GWs. 



2.2 Gravitational waves from turbulence and bubble collisions 

Let us first consider bubble collisions. A simulation with a large number of bubbles was 
carried out in Ref. [27] using the envelope approximation. This approximation neglects 
the overlap regions of colliding bubbles and follows only the evolution of the uncollided 
bubble walls (assuming thin fluid profiles). In Ref. [27] the bubbles were assumed to 
nucleate with a rate per volume and time given by 

r(t) = r(T i )exp[p(t-t i )], (8) 

and to expand with a constant velocity v w . The result for the peak frequency and intensity 
from the simulation is [27] 

f * ~ L6Xl ° ^1.8-0.1^ + ^ \m lOOGeVtf? (9) 
,coii nM |100V /3 0.11*4 ( PK \ 2 fH^ 2 



"r = o, 3i -j ^ Xfl[j) « R . (10) 

This result agrees with Eqs. (6) and (7), except for a slight difference in the dependence 
on v w . This can be seen by assuming L$ ~ 2v w j3~ 1 , since the time scale in this simulation 
is given by (3. The parameter /3 does not depend on details of the dynamics of the phase 
transition and is relatively easy to estimate for a given model. According to Eq. (8), we 
have 

/3 = r/r. (ii) 

The temperature decrease rate is governed by the Hubble rate, dT/dt ~ —HT. Hence, 
we have 

l = -™. (12) 

H TdT K J 

Concerning turbulence, one expects that eddies of a given scale L s will generate GWs 
with frequency given by /* ~ 1/As an d energy density given by Eq. (7). The size 
distribution of the eddies, as well as the energy distribution of the turbulence, is difficult 
to determine. In general, a single stirring scale Ls is assumed in the calculation. Below 
this scale, a Kolmogorov spectrum is established, according to which eddies of a given 



1 The ratio Pk/p is related to the generally used parameters k (the efficiency factor) and a (defined 
in section 3.1) by (pk/p)* = «a/(l + a) [notice that (pk/pr)* — H- 
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size break into smaller ones. This generates a cascade of energy which ends at a much 
smaller scale, related to the viscosity of the fluid. Above the stirring scale, the spectrum 
is determined by causality. In Ref . [13] , these two behaviors were assumed on each side of 
the stirring scale. In this case, the result for the peak agrees with the estimate (7). More 
recently [32], turbulence was modeled using a smooth interpolation between the large 
scale behavior and the small scale one. Besides, the fact that turbulence lasts for several 
Hubble times was taken into account. According to these results, the peak frequency is 
shifted to / p * ~ 3.5/ L s . Using the bubble size, L s ~ 2R b , Eq. (6) gives 

fT h = 2.7 x 10" 5 Hz ( -^-) V6 ( — ^ — L_ . (13) 

Jp V100/ \100 GeV J H*Rb y ! 

A fit to the result of Ref. [32] for the GW spectrum is given in Ref. [44] (see also [9]). The 
parametric dependence changes with respect to the estimate (7). For the peak intensity 
we have 

^ b = o.63 ( i {RbH l* ) «*■ (14) 

Notice that, since (RbH)* < 1, the dependence on the spatial scale in Eq. (14) is practi- 
cally that of Eq. (7), ^ urb oc (R b H)l 

2.3 Characteristic size scales 

As pointed out in Ref. [20], deciding which bubble size is relevant for turbulence is a 
major source of uncertainty. Bubbles of different sizes are present at the collision time. 
The bubbles of a given size stir up the fluid at that size scale, hence producing eddies of 
that scale. Larger bubbles in principle generate larger eddies and a larger GW intensity, 
whereas smaller bubbles in principle generate smaller eddies, but are more abundant. In 
any case, it is not clear which will be the turbulence spectrum in the case of several stirring 
scales. Although last nucleated bubbles are more abundant, they act during a shorter time 
and, furthermore, are probably "eaten" by larger bubbles. In our numerical calculations 
we shall use the size of the largest bubbles at percolation. Notice that, in contrast, there 
is no ambiguity in the production of GWs through bubble collisions, since the numerical 
fit (9)-(10) was obtained as a function of the parameters used in the simulation [27]. 
There is only some arbitrariness in the temperature at which the parameter j3 should be 
calculated, since in a real phase transition /3 is not a constant. 

We can obtain a simple estimate of the bubble size dispersion if we assume a nucleation 
rate of the form (8) and a constant wall velocity, and we neglect bubble overlapping. The 
bubbles which nucleated at time t have a radius Rb(t,t p ) = v w (t p — t) at the percolation 
time tp, and occupy a volume dV oc T(t)Rb(t,t p ) 3 . Thus, for the largest bubbles, which 
nucleated at time ij, we have 

Rma,x ~ v w(tp ~ ti)- (15) 

On the other hand, the bubbles which occupy the largest volume are given by the condition 
d(YR\)/dt = 0. Using Eq. (11), we see that the size Ry corresponding to the maximum 
of the volume distribution is given by 

R v = 3v w //3(t v ). (16) 



6 



This equation can be solved, taking into account that the nucleation time ty is related 
to Ry through Ry = v w (t p — ty). The time ty thus estimated should be between the 
nucleation time of the first bubbles U and the percolation time t p . On the other hand, 
the duration of the phase transition is usually assumed to be At ~ which gives for 
the largest bubbles 

(17) 

This approximation thus gives -R max ~ Rv- However, one expects Ry <C R ma , x due to the 
rapid variation of the nucleation rate. One could use the value (3~ 1 (t p ), which is larger 
than /3~ 1 (ty), in Eq. (17). However, this will not solve the problem. As we shall see in 
the next section, the parameter j3 does not have a huge variation between tj and t p [this 
validates the approximation (8)], except for very strong phase transitions. In general, the 
variation of /3 will not even compensate the factor 3 between Eqs. (16) and (17). For very 
strong phase transitions, on the other hand, /3(t p ) may become negative. This indicates 
that the widely used approximation At ps for the duration of the phase transition 
should be refined for this application. 

The estimate At ~ few//3 was obtained in Ref. [16], where in fact few = log(M/m), 
with M 3> 1 and m< 1. Notice that such a numerical factor in the bubble radius may be 
important. Indeed, the peak frequency (13) is proportional to l/Rb, whereas the intensity 
(14) is approximately 2 proportional to i? 2 . We can estimate the factor as follows. The 
initial nucleation time is that at which there is a bubble in a Hubble volume. Roughly, 

H~ 3 r F(t)dt ~ 1. (18) 



For a nucleation rate of the form (8) we obtain 3 F{Ti) ~ H 3 /3. The percolation time is 
roughly given by the condition 

ftp 47]- 

J -jv 3 w (t p -tfr(t)dt~i, (19) 

which yields 8ttvI F(ty exp[/3(t p - U)} ~ (3 4 . Since F(fy ~ H 3 (3, we have 

t p -t t ~3\og(J^j x/T 1 . (20) 

As we shall see, Eq. (20) is indeed a good approximation. Hence, we see that i^ max / Ry ~ 
\og(/3/H). If the log is of order 1, then we have t p — U ~ However, this is not in 
general the case. In general we have j3 3> H, since the bounce action varies very quickly. 
The value of (3/H is usually assumed 4 to be j3/H ~ 100, which gives t p — tj > 10/(3. 
Indeed, as we shall see in the next section, in general we have R max /Ry > 10. If the 
largest bubbles are relevant for turbulence instead of those of maximum volume, then the 
approximation Rt> ~ v w f3~ 1 in Eq. (14) leads to an underestimation of the GW signal 
from turbulence by at least two orders of magnitude. 



2 For R b H < 1, the deviation from the law fi* urb oc R% is at most a 2%. 
3 We see that the usual rough estimate r(Tj) <~ H A is valid only if j3 <~ H . 
4 As we shall see, /3 can depart significantly from this value. 
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3 Phase transition dynamics 



As bubbles expand, latent heat is released at the phase boundary. Part of this energy 
raises the temperature of the plasma, and another part is converted into kinetic energy 
in bulk motions of the fluid. The system we consider consists of the fluid and the Higgs 
field 0. All the thermodynamic quantities (energy density, pressure, etc.) are derived 
from the free energy J-"(0, T). In a range of temperatures around the electroweak scale 
T ~ 100 GeV, the high-temperature minimum of J 7 , = 0, coexists with a symmetry- 
breaking minimum m (T). The minima are separated by a barrier. In the unbroken- 
symmetry phase, the free energy density is given by J-+(T) = J-~(0,T), whereas in the 
broken-symmetry phase, it is given by F-(T) = J r (0 m (T),T). For a given temperature 
T, the pressure in each phase is given by p±{T) = —J-±(T), and the energy density is 
given by p±(T) = J-±(T) — TdJ r ±(T) / dT . The critical temperature is that for which 
•7-+(T c ) = F_(T C ), and the latent heat is defined as L = p + (T c ) — p_ (T c ). 

3.1 Bubble wall velocity and fluid profiles 

For hydrodynamic considerations we can assume an infinitely thin wall (see, e.g., [35]), 
such that the temperature and the fluid velocity are discontinuous at the interface. Con- 
sider a stationary wall which is locally moving in the x direction. We call T + and T_ 
the temperatures just in front and just behind the wall, respectively. The continuity 
conditions for energy and momentum fluxes give the relations [45] 



where v± are the values of the fluid velocity v on each side of the wall, in the rest frame 
of the wall, 7 = l/y/l — v 2 , w = p + p is the enthalpy density, and we use the notation 
p + = p + (T + ), p_ = p_ (T-), etc. These equations give v + as a function of v— The 
solutions have two branches, called detonations and deflagrations. For detonations the 
incoming flow is faster than the outgoing flow (|i>+| > |f-|)- The value of \v + \ is supersonic 
in all the range < |i>_| < 1, and has a minimum at the Jouguet point \v_\ = c s , where 
the speed of sound c s is given by c 2 s (T) = dp /dp. The minimum value of \v + \ is called 
the Jouguet velocity v j et . For deflagrations we have |i>+| < |i>_|, and |i>+| has a maximum 
value v j ef < c s at the Jouguet point |i>_| = c s . The hydrodynamical process is called weak 
if the velocities v + and t>_ are either both supersonic or both subsonic. Otherwise, the 
hydrodynamical process is called strong. 

There can also be discontinuities away from the phase-transition interface, which are 
called shock fronts. In this case Eqs. (21) and (22) still apply, only the equation of state 
relating the variables w, p, and T is the same on both sides of the discontinuity. As a 
consequence, the solution is simpler. The shock front is always supersonic. 

A macroscopic equation for the friction force exerted by the plasma on the bubble wall 
is usually obtained by introducing a phenomenological damping term and then integrating 



w + ^\v 2 + +p + = w-^tvt+p-. 



(21) 
(22) 
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the equation of motion for the Higgs field. One obtains 5 [36] 

P+-P--\ (s + + s_) (T+ - T_) + | (K| 7+ + |u_| 7 _) = 0, (23) 

where s = w/T is the entropy density and rj is the friction coefficient, which can be 
obtained from a microphysics calculation as explained in section 4. The various ther- 
modynamical variables are not independent, so Eqs. (21), (22) and (23) have only four 
unknowns, namely, the velocities v± and the temperatures T±. Besides, the temperature 
T + can be determined from the temperature T Q outside the bubble, which is known from 
the dynamics of the phase transition (see below). 

Out of the phase transition front, the fluid velocity profile (in the reference frame of 
the bubble center) depends on the symmetry of the bubble. This issue was discussed in 
Ref. [40]. The total amount of energy transmitted to the plasma is qualitatively and 
quantitatively similar for different wall geometries. We shall consider planar walls, which 
are simpler and allow to obtain analytical expressions (notice that, as bubbles collide, 
any previous symmetry is lost). For a stationary wall moving with velocity v w , there 
is no characteristic distance scale in the fluid equations. As a consequence, it is usual 
to assume the similarity condition [45] , namely, that the temperature and velocity of the 
fluid depend only on £ = x/t. For the planar symmetry case, we have for the fluid velocity 
(see e.g. [40]) 

v 1 = 0, (24) 

where a prime indicates derivative with respect to £. This equation gives either constant 
solutions v (£) = const, or a "rarefaction wave" solution 

Similarly, the enthalpy profile is given by the equation 

w ' f 1 , i~i 2 , 

— = l ~2 + 1 J ~\ — r 1 v ' 

w \ c i ) 1 -i v 

which can be readily integrated if v is a constant or the rarefaction solution (25). The 
fluid velocity and temperature profiles are thus constructed with these solutions, using 
the matching conditions (21) and (22) and appropriate boundary conditions. 

The usual boundary conditions consist of a vanishing fluid velocity far behind the 
moving wall (at the center of the bubble) and far in front of the wall, where information 
on the bubble has not arrived yet. We shall refer to the temperature far in front of the wall 
as the "outside" temperature T . The initial value of T a is the temperature T n at which 
the bubble nucleated, but T Q will change due to the adiabatic expansion of the universe 
or the presence of other bubbles. To be consistent with the similarity condition, however, 
T Q should be a constant T Q = T n , so that the wall velocity would also be a constant. We 

5 See [42] for a discussion on the validity of this equation. 
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(25) 



(26) 



shall assume that T Q changes slowly enough to allow the wall to be always in stationary 
motion. 

Three kinds of solutions for the wall velocity and fluid profiles are possible (for a recent 
discussion see [42]), namely, a weak detonation, a "traditional" weak deflagration, and a 
supersonic Jouguet deflagration. Let us denote v± the fluid velocities just in front and 
behind the wall, i.e., 

_ _ v± + v w 
1 + v w v± 

The wall is at the position £ w — v w . For the detonation solution, the wall is supersonic 
and the fluid in front of it is unperturbed. Therefore, the fluid velocity v + vanishes and 
we have v w = \v+\. It turns out that the detonation solution can only be weak or, as a 
limiting case, Jouguet. Behind the wall, the fluid velocity is a constant v = V- up to a 
point £o which lies between c s and £ w . At £ = £o the fluid velocity matches the rarefaction 
solution (25). Continuity implies 

6 = (28) 

1 + V_C S 

The rarefaction solution vanishes at £ = c s , and we have v = for £ < c s . For the 
traditional deflagration solution, the fluid behind the wall is at rest, so v- = and 
v w = \v-\. Again, this solution can only be weak or, at most, Jouguet. Therefore, the 
wall is subsonic. The fluid velocity in front of the wall is a constant v = v + up to a shock 
front where the profile ends, at a point £ s h > c s determined by the shock discontinuity 
conditions (see below). Beyond the shock, the fluid is still unperturbed. Finally, the 
supersonic deflagration is a Jouguet deflagration. In this case, the condition v- = of 
the traditional deflagration is replaced by the Jouguet condition t>_ = — c s , and we have 

v- = (v w - c a )/ (1 - v w c s ). (29) 

The wall velocity is always supersonic and the fluid velocity behind the wall is given by 
the rarefaction solution (25) between c s and £ w . In front of the wall the fluid velocity is a 
constant v — v + between £ w and £ s h- In the limit £ w — c s , there is no rarefaction wave and 
the solution matches the traditional deflagration. As £ w increases, the shock front and 
the phase-transition front become closer. As £ m reaches the Jouguet detonation velocity 
Vj et , the shock wave disappears and the solution matches the detonation. 

Equations (21), (22), and (23) for the wall velocity, and Eqs. (25) and (26) for the 
profiles, can be solved once the equation of state (EOS) of the system is known. It is 
convenient to approximate the model by using the bag EOS 

jr + (T) = -a + T 4 /3 + e, T_ (T) = -a_T 4 /3, (30) 

which corresponds to having only radiation and vacuum energy in the symmetric phase, 
and only radiation in the broken-symmetry phase. This simplification allows to find 
analytical expressions for the solutions. In this model the latent heat is given by L = 4e 
and the speed of sound is a constant, c s = 1/V3. It is customary to express the results 
as functions of the variable a = ej (a + T 4 ) (which gives the ratio of the vacuum energy 
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density to the energy density of radiation). As discussed in Ref. [37], for applications it 
is convenient to use the latent heat L instead of e. Therefore, we define the parameters 

L L L 

4a+T c 4 ' + 4a + T 4 ' 4a + T Q 4 ' v ; 

corresponding to the critical temperature T c , the temperature just in front of the bubble 
wall T + , and the temperature far in front of the wall T . The solutions for the wall 
velocity and fluid profiles will depend only on a c , a Q , and r]/L. The temperature T Q is 
the boundary condition for the temperature profile. The corresponding enthalpy density 
is given by 

w = \a + T* = A (32) 



The matching conditions relating the values of W-, w + and w are given by Eq. (21) at 
the phase-transition discontinuity and the equivalent equation for the shock discontinuity. 

Using the equation of state (30) in Eqs. (21) and (22) we obtain the relation between 
v + and V-, which for this model depends only on the parameter a + [46], 



v + = — . (33) 

The plus sign corresponds to detonations and the minus sign to deflagrations. The friction 
equation can also be expressed in terms of t> + , and a + , 



with 



1 - 3v + v^ 3 V s+/\ T. 



s + a + \T + ) T + 



a + ( ' 1 + v + v. 

— 1 - OLa 



1/4 

(35) 



a_ \ 1/3 — v + V- 

The ratio a_/a + is given by a_/a + = 1 — 3a c . From Eqs. (33)-(35) we can find the 
velocities v + and f_ as functions a + . The relation between a + and a a depends on the 
type of hydrodynamic solution. For detonations, the wall velocity is given by v w = 
—v + and the temperature T + is just the outside temperature, hence a + = a Q . For 
deflagrations, the temperature T + is related to T Q through the matching conditions at the 
shock discontinuity, which for the bag EOS are given by 

1 Vl 3T Q 4 + T\ 

ViV2 = s> "WTTr (36) 

where V\ is the velocity of the outgoing flow in the reference frame of the shock, and v 2 
that of the incoming flow. In the frame of the bubble center, the fluid velocity in front of 
the shock vanishes. Hence, the shock velocity is given by i> s h = — 1>2, and we obtain 




^sh = t + \ T +o> ( 37 ) 
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where Vi is the fluid velocity behind the shock. In the shock-wave region the fluid velocity- 
is a constant. As a consequence, we have v\ — v+, which gives 

v+ - Vw = VS(ao - a+) , . 

1 - v+v w v /(3a + a + ){3a + + a )' 

For traditional deflagrations, we have i>_ = —v w , so Eq. (38), together with Eq. (33), 
can be used to obtain function of a a . For Jouguet deflagrations, V- = —l/y/3 is 

fixed, so Eq. (33) alone gives v + as a function of a+, i.e., v + = Vj ci (a + ). In this case, 
Eq. (38) gives the wall velocity as a function of a + and a Q) and Eq. (34) can be used to 
eliminate a + . 

There is in general a solution for any set of parameters. As a matter of fact, there can 
be more than one. In such a case, only one of them will be realized in the phase transition. 
This issue is discussed in Ref. [42], where the semi-analytical solutions of Eqs. (33)- (38) 
are compared with a numerical calculation [47]. As a general rule, the weak detonation 
is more stable than the Jouguet deflagration, and the latter is more stable than the 
traditional deflagration. As an example, we show in Fig. 1 the solutions that are realized 
as final stationary states, as a function of the friction, for the values of the parameters 
considered in Ref. [42] . For large values of the friction we have weak deflagrations. When 

1.0 
0.8 
0.6 
0.4 
0.2 




°O01 0.05 0.10 0.50 1.00 5.00 10.00 

TJ/L 

Figure 1: The wall velocity (solid line) and shock velocity (dashed line) as functions of 
the friction, for a c = 4.45 x 1CT 3 and a Q = 7.06 x 10~ 3 . 

the speed of sound is reached, i.e., at the Jouguet point, the traditional deflagration 
matches the supersonic (Jouguet) deflagration (notice the discontinuity in the derivative 
of the curve, which is due to the change of hydrodynamical solution). For a lower value of 
the friction, the detonation solution appears. Since this solution is the stable one, there 
is a jump in the wall velocity as a function of r\. The dashed line indicates the velocity of 
the shock front. We have v s h ~ c s for subsonic deflagrations. 
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3.2 Energy injected into the plasma 

The energy released at the phase transition (i.e., the latent heat) reheats the plasma and 
causes bulk motions of the fluid. The generation of GWs requires the spherical symmetry 
to be lost. This happens once bubble walls or shock fronts collide. So far we have 
considered fluid profiles for stationarily moving walls, since it is hardly possible to know 
the profiles during bubble collisions. This is irrelevant for the envelope approximation, 
which only takes into account the motion of uncollided walls (and assumes thin profiles). 
To estimate the average energy which goes into the formation of turbulence, we may 
calculate the energy density of the fluid just before the profiles meet, and assume that, 
after the fronts collide, this energy gets redistributed throughout the space occupied by 
the bubbles. With this assumption, we only need to consider the average energy density 
for isolated bubbles. Notice that, although there are bubbles of different sizes (because 
they nucleated at different times), the wall velocities and fluid profiles depend only on 
the temperature T (which is the same for all bubbles). 

The kinetic energy density of the fluid is given by pkm — wv 2r y 2 . Let us first consider the 
subsonic deflagration. For planar walls, the kinetic energy density is a constant between 
the wall and the shock front, pkin = ^+£+7+, and vanishes elsewhere. Therefore, the total 
energy is proportional to R s h—Rb, where Rb and i? s h are the positions of the bubble and the 
shock front, respectively. To calculate these positions we should integrate the respective 
velocities. However, the profiles were calculated using the similarity condition and, for 
consistency, we must consider again this approximation 6 . Thus, we have Rb = £ w Atb, 
Rsh = £ s hAib, where At b is the time during which the wall has been moving, £ w is the 
wall velocity calculated from Eqs. (33)- (38) and £ sh is the shock front velocity given 
by Eq. (37). Assuming that, after the shock fronts meet, the energy which was initially 
concentrated in front of the bubble wall gets distributed in the whole volume proportional 
to £ s hAt b , the average kinetic energy density in the fluid is given by 

p K = w+v 2 + ^ 2 + ^— — — (subsonic deflagrations), (39) 

which is the same for bubbles nucleated at different times. 

For detonations, the kinetic energy density is concentrated between c s and £ w . Between 
c s and £ there is the rarefaction wave, and between £ and £ w the fluid profiles are 
constant. The integration of the kinetic energy density in the rarefaction region was done 
analytically in Ref. [40]. In the case of detonations, bubble collisions and turbulence 
begin when the bubble walls meet (since there are not shock fronts). Dividing the total 
kinetic energy by the volume of the bubble, we have 



PK = W. 

where 



, 3/ /-\3r fl-v\vSf(C w )-f(c s ) 



~2 ~2 Su> ~~ SO . <J f n AT\ 



1 + v. 



(detonations), 
(40) 



i -<ey IV3 



y/3 2 ' 



(41) 



6 Our numerical calculation shows that the velocity generally varies by at most a 30% before colliding. 
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and 2 Fi is the hypergeometric function. 

The profile of a supersonic deflagration consists of a shock wave in front of the wall 
and a rarefaction wave behind it. The average kinetic energy density is thus given by 

3 f l~Cw \^ f(Cw) -f(Cs) , - 2 ~2 U - ( in i- \ 

Pk —w- I — 1 h%w + 7 + — (supersonic deflagrations). 

4 V 1 + <Lw J ?sh £sh 

(42) 

3.3 Bubble nucleation, expansion and percolation 

In principle, as soon as the temperature descends below T c , bubbles begin to nucleate 
and expand. However, at the beginning there will be too few bubbles. The "onset" of 
nucleation is usually defined as the time at which there is already one bubble in a Hubble 
volume. We shall take this as the nucleation time of the "first" bubbles. These will be 
the largest bubbles in the system, thus setting the characteristic wavelength of the GWs. 

On the other hand, bubble collisions can in principle begin as soon as there is a non- 
vanishing probability of having a couple of bubbles in a causal volume. However, at the 
beginning the bubbles will be too few and too small, hence their collisions will be very 
unlikely. Bubbles will effectively begin to meet and collide once their density and size have 
become large enough. At first, there will form clusters of a few bubbles, and then larger 
and larger clusters. Percolation occurs when a cluster of infinite size spreads through the 
medium (equivalently, when there is a cluster spreading from side to side in a large box, 
say, of Hubble size). Percolation has been studied numerically for spheres (of equal size) 
in a large box. With the spheres distributed at random and allowing overlapping, an 
infinite chain is established when the fraction of space covered by spheres is 0.29 [48]. We 
shall assume that, at this moment, most bubbles are already colliding. 

The nucleation of bubbles [49, 50] is governed by the three-dimensional instanton 
action 



oo 

2, 



S3 = 47T / ;-(// 

where 







1 " - 2 



- 2 m + V T Mr)) 



(43) 



Vt(0)^^(0,T)-^(O,T). (44) 

The bounce solution of this action, which is obtained by extremizing S3, gives the radial 
configuration of the nucleated bubble, assumed to be spherically symmetric. The action 
of the bounce coincides with the free energy of a critical bubble in unstable equilibrium 
between expansion and contraction. The solution obeys the equation 

d 2 (j) 2 d(j) _ dV T 
dr 2 r dr d<p ' 

with boundary conditions 

^(0) = 0, lim <j>(r) = 0. (46) 

dr r^oo 

The thermal tunneling probability for bubble nucleation per unit volume per unit time is 
[50] 

r(T) ~ A(T) e" 5a(T)/T , (47) 
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with A(T) = [S 3 (T)/(2irT)] 3/2 T 4 . At the critical temperature, S 3 diverges and the nucle- 
ation rate vanishes, whereas at the temperature at which the barrier between the minima 
of J 7 disappears, S3 vanishes and the nucleation rate becomes extremely high, T ~ T 4 . 
We define the nucleation time tj of the first bubbles by the condition 

V H n(U) = 1, (48) 

where Vh = H~ 3 is the Hubble volume, and the number density of bubbles is given by 

'ait') 1 3 



n{t) = I dt'T(T(t')) 

Jtr. 



a(t) 



(49) 



where t c is the time at which the Universe reached the critical temperature T c . The scale 
factors take into account the fact that bubbles which nucleated at time t' with a number 
density dt'T (T(t')), get diluted until time t due to the expansion of the Universe. At this 
initial stage, the temperature variation is determined by the adiabatic expansion equation 
dp + = —3w + da/a. Hence, we have 

dF+/dT da 

dT = - 3 d^W 2 ^ m 

The evolution of the scale factor is given by the Friedmann equation 

— = Hdt, (51) 
a 

with the expansion rate given by 

H = y/8irGp + /3. (52) 

If the high-temperature phase is comprised only of radiation and vacuum energy, then Eq. 
(50) becomes dT = —Tda/a, which gives the well known result dT/dt = —HT. However, 
some of the models we consider in the next section have particles with masses ~ T in the 
+ phase. 

When bubbles begin to nucleate, the energy density is no longer homogeneous. On 
the one hand, even if the temperature were homogeneous, the internal energy of the — 
phase is lower than that of the + phase. On the other hand, temperature gradients 
arise due to the latent heat released at the interfaces. The Hubble rate is thus governed 
by the average energy density. In fact, energy conservation implies that the released 
energy compensates the decrease in the broken-symmetry phase. As a consequence, the 
average energy density will not depart significantly from p + (T ), which decreases due 
to the adiabatic expansion. Therefore, we shall use Eq. (52) still in the presence of 
bubbles. This is a good approximation for the stages of the phase transition we are 
interested in (i.e., up to the percolation time). In any case, the scale factor will not 
change significantly during the inhomogeneous stage. We have checked numerically that, 
as soon as the broken-symmetry regions become barely appreciable, say, at a time t when 
the fraction of volume occupied by bubbles reaches a value fb ~ 1CT 2 , the phase transition 
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is already happening so quickly (due to the extreme behavior of the nucleation rate) that 
the percolation fraction fb ~ 0.3 is achieved in a time St <C i — t c . 

Assuming a homogeneous nucleation throughout the symmetric-phase regions, the 
fraction of volume occupied by bubbles is fa — 1 — f u , where f u is the fraction of space in 
the unbroken-symmetry phase, given by [51] 



The radius of a bubble that nucleated at time t' and expanded until time t is given by 



where R n is the initial radius of the nucleated bubble, which immediately becomes negligi- 
ble in comparison to the second term in Eq. (54). We have used the notation T' Q = T (t'), 
a' = a(t'), etc. Notice that, at a given time t, all bubble walls move with velocity v w (T (t)), 
where T Q (t) evolves according to Eq. (50). We shall assume that, as the temperature de- 
creases, the hydrodynamics instantaneously adjusts to a stationary solution. Moreover, if 
the stable type of stationary solution changes, e.g., from a deflagration to a detonation, 
we approximate the velocity variation by a jump. The factors of a in Eqs. (53) and (54) 
take into account the fact that the number density of nucleated bubbles gets diluted and 
the radius of a bubble gets stretched due to the expansion of the Universe from t' to t. The 
exponent in Eq. (53) would give a naive result for f b assuming a homogeneous nucleation 
rate throughout space (including the broken-symmetry regions). Thus, Eq. (53) avoids 
overcounting of overlapping or nested bubbles. This result is obtained by considering the 
probability that a given point in space lies outside of any bubble (this is why a bubble 
that nucleated in the broken-symmetry region does not contribute to fb even though it 
contributes to the exponent). 

However, Eq. (53) assumes that the nucleation is homogeneous in the symmetric 
phase, with a rate T(T ). In fact, temperature profiles may cause considerable inhomo- 
geneities in the nucleation rate. Consider bubbles which have not yet interacted with 
each other. If the hydrodynamic solution is a detonation, then the temperature in the 
symmetric phase is just T = T + = T G , and Eq. (53) does indeed hold. On the other hand, 
in the case of deflagrations there is a reheated zone in front of the bubble walls (T + > T Q ). 
Since the nucleation rate is extremely sensitive to temperature, bubble nucleation is ef- 
fectively turned off in such regions. Therefore, we can assume that the nucleation rate 
vanishes, not only inside the bubbles, but also in the shock-wave regions, whereas it is 
given by T (T Q ) beyond the shock fronts. Equation (53) does not take into account this 
fact, and must be modified in order to avoid bubble overcounting. We accomplish this by 
considering, instead of the fraction of volume / s h occupied by "shock-front bubbles" , 
which is obtained by replacing the bubble radius Rb in Eq. (53) by the shock 
front radius R sh calculated from the shock front velocity i> s h instead of v w . Moreover, 
in the deflagration case we are not interested in fb but in / s h, since turbulence begins as 
soon as the shocks collide. 

We shall follow the evolution of the phase transition up to the percolation time t p , 
which we define as the moment at which the fraction of volume occupied by bubbles (in 




(53) 




(54) 
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the case of detonations) or by shock bubbles (in the case of deflagrations) reaches the 
value 0.3. We will solve Eq. (45) iteratively by the overshoot-undershoot method, and 
we will integrate Eq. (43), Eq. (49), and the set of Eqs. (50)-(54) numerically (see Ref. 
[26] for details). The relevant quantities, such as the temperature T (t p ) or the size of the 
bubbles (which roughly goes with t p — ti) are not much sensitive to the definitions of ti 
and t p . Indeed, these definitions involve the bubble number density n and the fractions of 
volume fb or / s h which, due to the extreme behavior of T with t, change by many orders 
of magnitude in the characteristic time t p — t^. As a consequence, changing the values of 
n or / that one uses to define ti and t p (even by an order of magnitude) introduce very 
small differences At itP <C t itP . We have checked this issue numerically. 

In order to compute the GW signal from bubble collisions, we only need to evaluate 
the wall velocity, the kinetic energy of the fluid, and the parameter j3 defined in Eq. (11). 
Since the dynamics of nucleation is dominated by the variation of S3, we can neglect the 
variation of the prefactor in Eq. (47). Therefore, we have 

This parameter must be computed at some characteristic temperature. It is not clear 
whether this temperature should be chosen close to Tj or rather to the later temperature 
T p . In the simulation which gives the fit (9)-(10), (3 is a constant [27]. In Fig. 2 we show 
the values of (5 both at the initial time and at the percolation time for one of the models 
considered in the next section. We see that the difference is a factor of 0(1) (generally 
less than 2), except for extremely strong phase transitions (which are those very close to 
the maximum value of the parameter h s in Fig. 2). In general, this difference 7 is quite 
smaller than other uncertainties (see below), and we shall use the value /3(Tj). 

In Fig. 2 we also show the value of S 3 /T at 7* and T p . We see that S 3 (Ti)/Ti takes 
values around the well known estimation S3/T ~ 140, which can be deduced from Eqs. 
(18)-(20). Notice also that the value of /3 departs in general from the usual assumption 
j3 ~ 100H . This assumption is based on the argument that the time scale for change in 
the nucleation action should be comparable to that in which the temperature changes [16] , 
which gives f3/H ~ S 3 /T. We see that this estimate is too crude. We also see that the 
approximation At ps 3 log(/?/iJ)/3~ 1 gives a very good estimate of the time At = t p — U. 
In contrast, the estimation At ~ is at least an order of magnitude smaller than At. 

In order to compute the contribution of turbulence to the generation of gravitational 
waves, we shall calculate the average energy density in bulk motions of the plasma at the 
percolation time. We also need to know the typical size of the bubbles. As discussed in 
section 2, there is an arbitrariness in the appropriate size scale L s , and we shall consider 
the largest bubbles. We shall calculate the radius Rb(ti,t p ) of the largest bubbles in the 
case of detonations, or the radius R s h(ti, t p ) of the largest shock-front bubbles in the case of 
deflagrations, by integrating the corresponding velocity v w or v s h from ti to t p . Regarding 
the widely used approximation Rf, ps v w f3~ 1 for the largest bubbles, as we have seen, 
R b ^ v w 3 \og{(5 / H)^ 1 will give a much better approximation. In Fig. 3 we show the 

7 The difference between /3(Tj) and /3(T P ) becomes larger near the maximum value of h s . However, 
for very strong phase transitions it is convenient to evaluate j3 at Ti since this parameter may become 
negative for small temperatures (see, e.g., Fig. 5 of Rcf. [22]). 
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Figure 2: The inverse time scales and the bounce action, for an extension of the SM with 
a complex scalar singlet with coupling h s to the Higgs and invariant mass fi s = 0. Red 
lines are calculated at t — ti, black lines at t — t p , and the blue line corresponds to t p — ^. 

largest bubble radius at the percolation time, together with the two estimations. For this 
particular example, the approximation At = underestimates the radius by a factor 
30 in most of the parameter range considered in the figure. Hence, this approximation 
will underestimate the amplitude of the waves by a factor « 900. Figure 3 also shows 
that considering a constant velocity v w = v w (Ti) is in general a good approximation, as 
well as (3 = /3(T;). 

4 The electroweak phase transition 

In the SM, the electroweak phase transition is only a smooth crossover [33]. However, 
many extensions of the model give a first-order phase transition. For simplicity we shall 
consider models with a single Higgs field, or models for which a single Higgs provides a 
good approximation. 

4.1 Effective potential and bag parameters 

Our theory will consist of a tree-level potential 

V (0) = -m 2 2 + (56) 
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Figure 3: The radius of the largest bubbles at t = t p and the different approximations, 
for the same model and parameters of Fig. 2. Black lines correspond to the bubble wall 
and red lines to the shock wall. The parameter /3 is calculated at the initial temperature 
T 
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for the background Higgs field 0, defined by (H°) = <f>/y/2. The vacuum expectation 
value (vev) of is given by v = ^/2/\m = 246 GeV, and A fixes the Higgs mass, m 2 H = 
2Xv 2 . Imposing the renormalization conditions that the minimum of the potential and 
the mass of do not change with respect to their tree-level values [52], the one- loop 
zero-temperature effective potential is given by V ((f)) = V (0) + V 1 ((f)), with 



+ c, (57) 



where the upper and lower signs correspond to bosons and fermions, respectively, gi is the 
number of d.o.f. of the particle species i, ((f)) is the 0-dependent mass, and the constant 
c is chosen so that the energy density vanishes in the true vacuum at zero temperature 
[37], V (v) + Vi(v) = 0. The free energy (finite-temperature effective potential) to one- loop 
order, including the resummed daisy diagrams, is given by 

T{4>,T) = V ((f)) + V x ((f)) + J\(0,T), (58) 

where the finite-temperature corrections are given by [53] 

J\(0,T) = ^±|^jT^ 2 log l T ^{-\J^ + mj (0)/T^ 

+ E K (0 " Ml ((f))] . (59) 

bosons 

Here, Aii is given by Aif ((f)) = mf ((f)) + Hj (T), where Hj (T) is the thermal mass. 
The last term receives contributions from all the bosonic species except the transverse 
polarizations of the gauge bosons. Some of the masses in the gauge sector are gauge 
dependent. We use the Landau gauge and, thus, we consider only the transverse and 
longitudinal polarizations of the gauge bosons 8 . 

We shall consider in general Higgs-dependent masses of the form 

m 2 ((f)) = h 2 (f) 2 + fil (60) 

For Hi T, the contribution of the species i to the energy density of the unbroken- 
symmetry phase is that of radiation, i.e., proportional to ^T 4 . Since this is true for most 
species, we have in general p + ~ -K 2 g if T A /30 + p vac , where g* is the number of relativistic 
d.o.f., and p vac = Vo(0) + Vi(0) is the false vacuum energy density. In such a case, the bag 
parameters are given by a + = 7i 2 g*/ '30 and e = p vac . In the general case, we define the 
thermal energy density by = p + — p vac and compute the bag parameters a c , a Q , etc., 
by 

a (T) = — r— -, (61) 

where the latent heat is given by L = T c (dT-jdT — dj r + /dT) \t c (which does not fulfil in 
general the bag relation L = 4e). 



8 Although physical quantities such as the latent heat or the amount of supercooling should not exhibit 
any gauge dependence, an inconsistent truncation of the perturbative expansion can introduce a nontrivial 
gauge dependence [54] . This would be important in a model for which the strength of the phase transition 
relics on the gauge fields, which is not the case of the models considered here. 
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4.2 Friction force 



The friction was calculated in several studies of the microphysics [55, 56]. Some general 
approximations were derived in Refs. [25, 37]. The friction coefficient appearing in Eq. 
(23) receives contributions from thermal particles, i.e., those which obey the Boltzmann 
equation, and from infrared bosons, i.e., infrared excitations of bosonic fields, which must 
be treated classically. For masses of the form (60) we have, for thermal particles, 



E T]T T ^ (^/ T )] 2 W T ) 2 V^ d<P, (62) 



\ ^ :/:'' '' 

Vth 



where Vt is defined in Eq. (44), the limits of integration correspond to the minima of the 
free energy at T = T c , the function c\ is given by 



C1{X) = ^J X dV ^^^W (63) 

and f is an average interaction rate arising from the collision term of the Boltzmann 
equation. For the electroweak phase transition, f is typically ~ 1CT 2 . For infrared 
bosons, we have 

V, = £ 9j ^T f C b(m t /T) {<P/TY^V T d<P, (64) 

bosons < ^ >0 

where mo is the Debye mass, given by m 2 D = (ll/Q)g 2 T 2 for the W and Z bosons of the 
SM, and m 2 D = h 2 T 2 /3 for a scalar singlet. The integral in (64) has an infrared cut-off 
0o for small given by O = \/ L~ 2 — Jj^/h for /ij < L~ l , and O = for > L" 1 , 
where L w is the wall width. In the thin wall approximation, L w can be estimated as 
L w w d(j)/^W^. The function b is given by 

Each of these two contributions dominates in different parameter regions, and we can use 
V — Vth + Vir- Analytical approximations for i] th and i] ir in different limits can be found in 
Ref. [37]. 

The above expressions for the friction coefficient were derived in the small-velocity 
regime. Recently, the ultra-relativistic regime was considered in Ref. [39]. In this limit 
the friction does not depend on the velocity. The total force per unit area acting on a 
wall which is already propagating ultra-relativistically (with gamma factor 7 ~ 10 9 ) is 
given by 

F tot /A = p„-p + = -f_ + f + , (66) 

where J-{4>, T) is the-mean field effective potential. The latter is obtained by keeping only 
the quadratic terms in a Taylor expansion of the finite-temperature part of J- ((f), T) about 
= [39, 38], 



T{4>, T) = V (0) + V X (0) + $>. 2 (0) - m 2 (0)] ^ 



(67) 
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For the case of bosons with masses of the form m; = h^ (i.e., //j = 0, which gives stronger 
phase transitions), the last term in Eq. (67) becomes 



For fermions there is an additional factor 1/2. 

If the total force (66) is positive, then the wall can run away. This means that 
the bubble may undergo accelerated expansion instead of reaching one of the station- 
ary states considered in the previous section. As shown in Ref. [39], the bubble wall 
never runs away in a "fluctuation induced" first-order phase transition, i.e., a phase tran- 
sition which is first-order due to the thermal part of the potential. The typical example 
of a fluctuation-induced first-order phase transition occurs when the high-temperature 
expansion of T\{(j),T) has a cubic term — £V Tmf ((/>)/ 12-7T. This term is not present in 
the mean field potential. If the first-order character of the phase transition is due to this 
term alone, then in the mean field potential the broken symmetry minimum raises above 
the symmetric minimum. As a consequence, the force (66) is negative, which means that 
the wall cannot run away. An example of a model which may yield a strong enough 
first-order phase transition is a potential with tree- level cubic terms [39]. This is possible, 
e.g., in extensions of the Standard Model with singlet scalar fields. 

5 Numerical results 

The relevant SM contributions to the one-loop effective potential come from the Z and 
W bosons, the top quark, and the Higgs and Goldstone bosons. It is usual to ignore the 
Higgs sector in the one-loop radiative corrections. This should be a good approximation 
in extensions of the SM which include particles with strong couplings to 0. The <p- 
dependent masses of the weak gauge bosons and top quark are of the form hi(f>, with 
hi = rrii/v, where rrii are the physical masses at zero temperature. We shall ignore the 
longitudinal components of the weak gauge bosons, which are screened by plasma effects. 
Thus, the W and Z contribute corrections of the form (57), (59), with 4 and 2 bosonic 
d.o.f., respectively. The top contributes g t = 12 fermionic d.o.f. The rest of the SM 
particles have hi ^ 1 and only contribute a ^-independent term — it 2 g^g^T 4 / 90, with 
flight ~ 90. We shall consider several extensions of the SM, which may provide a strongly 
first-order phase transition. For all the models considered below, we use a Higgs mass 
rxiH = 125 GeV. The dependence of the relevant quantities on the strength of the phase 
transition is illustrated in Fig. 4. For further results on the wall velocity for these models, 
see Ref. [37]. For Results on the temperature and duration of the different stages of the 
phase transition, see Ref. [26]. In this paper we shall focus on the gravitational waves 
generated during the phase transition. 

5.1 Extra scalars 

The simplest extension of the SM consists of adding gauge singlet scalars [52, 57], which 
may range from a single field S [58, 59, 60] to several fields Si [61]. In general, these bosons 




(68) 
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constitute a hidden sector which couples only to the SM Higgs doublet through a term 
h 2 s H^H^2, Sf (assuming, for simplicity, universal couplings hi = h s ). The scalars may 
have £77(2) x U(l) invariant mass terms [i 2 S 2 and quartic terms A s £ 4 . For simplicity, 
we shall set A s = for our numerical calculations. We have checked that considering 
\ s 7^ does not introduce qualitative differences in the results. A negative value of n 2 
may enhance the strength of the phase transition. This fact is exploited in the case of 
the MSSM in the light-stop scenario, which we consider in the next subsection. Besides, 
adding real singlets allows cubic terms of the form (H^H) S or £ 3 , which cannot be 
constructed with Higgs doublets. The presence of cubic terms in the tree-level potential 
makes it easier to get a strongly first-order electroweak phase transition [58]. This may 
cause a significant increase in the wall velocity [37] and even the existence of runaway 
solutions [39]. As a consequence, tree- level effects may lead to an important GW signal 
from bubble collisions [9]. In order to study such tree- level effects, one should consider the 
full potential depending on the two fields H and £. This is out of the scope of the present 
paper, since our numerical calculations are based on a single- variable potential V ((/>). 
Therefore, we shall ignore the possibility that cubic terms exist in the tree-level potential. 
Thus, the contributions of the scalars to the free energy are of the form (57), (59), with 
m 2 ((f)) = h 2 <f) 2 + fi 2 and g s given by the number of real singlets. The thermal mass is 
given by IT = h 2 s T 2 /3 [59]. 

It is well known that the phase transition is more strongly first-order for larger numbers 
of bosons g s and stronger couplings h s (see Fig. 4). On the contrary, for high values of 
fi s the bosons decouple from the thermal plasma, and the phase transition becomes more 
weakly first-order. As a consequence, the wall velocity and the energy injected into 
macroscopic motions of the fluid generally grow with g s and h s and decrease with /i s . 
Indeed, increasing g s and h s , the separation between minima of the free energy increases, 
as well as the height of the barrier separating them. As a consequence, bubble nucleation 
effectively begins at lower temperatures Tj [22, 26]. For this particular model, the strength 
of the first-order phase transition has a strong dependence on the coupling h s . For large 
enough values of h s there is a barrier at zero temperature, and the nucleation temperature 
Ti is very small. As shown in Fig. 4, there is a value h s = /i max for which Tj falls to 
(upper left panel). Beyond this value the phase transition is too strong to overcome the 
supercooling stage and the Universe will eventually enter a period of inflation [22, 26]. 
This is reflected also in Figs. 2 and 3. The endpoint in the curves corresponds to /i max - 
Near this endpoint, the system remains stuck in the symmetric phase for a long time 
U — t c 3> H^ 1 (upper right panel of Fig. 4). Meanwhile, the temperature decreases to a 
value Tj <C T c . The same happens to the time At = t p — ti needed to arrive at percolation, 
as can be seen in Fig. 2. 

The generation of gravitational waves depends mainly on the wall velocity and the 
kinetic energy in bulk motions of the fluid. We show the values of these quantities at 
t — U in the lower panels of Fig. 4. The behavior at the percolation time is similar. 
We can distinguish a cusp in these curves, which indicates the passage from a value of 
h s for which the hydrodynamical solution is a weak deflagration, to a value of h s for 
which the wall moves as a Jouguet deflagration. Similarly, there is a jump indicating the 
passage from a Jouguet deflagration to a weak detonation. Notice that, although the wall 
velocity for the detonation is higher at the discontinuity point, the detonation is a weaker 



23 



0.0 

1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 




1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 



h 

S 



Figure 4: Several quantities calculated at the beginning of bubble nucleation, as functions 
of h s for g s = 2 and fi s = 0: the temperature and field mean value (upper left), the 
supercooling time ti — t c (upper right), the bubble wall and shock velocities (lower left), 
and the kinetic energy of the fluid (lower right). 
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hydrodynamical solution than the supersonic deflagration and causes a lower disturbance 
of the fluid. Therefore, the jump in the injected kinetic energy is negative. The local 
maximum at this discontinuity is due to the fact that the kinetic energy is maximal for 
Jouguet deflagrations [38, 40]. 

This behavior is reflected in the gravitational wave generation, as can be seen in Figs. 
5, 6 and 7. In Fig. 5 we plot the peak of the GW spectrum from turbulence, as a 
function of h s for different values of g s and /i s . As the parameters of the model are varied, 
the peak frequency and intensity change by several orders of magnitude. This variation 
includes frequencies in the sensitivity range of LISA and eLISA, / ~ 1 mHz. However, the 
peak sensitivity of LISA, li 2 fi G w ~ 10~ 12 , and that of eLISA, h 2 Qcw ~ 10~ 10 (marked 
with horizontal dotted lines in Fig. 5), are reached near the maximum values of h s . 
Unfortunately, for such high intensities the spectra do not peak at mHz frequencies. For 
values of h s which give mHz frequencies, the intensities are a few orders of magnitude 
below LISA's sensitivity. This is better seen in Fig. 6, where the peak value Q p is shown 
as a function of the peak frequency f p , together with the sensitivity curves 9 for LISA and 
eLISA. The curves for the predicted signal cross LISA's sensitivity curve at f p ~ 10~ 4 Hz. 
This happens very close to the endpoints. Reaching LISA's sensitivity thus requires to 
tune the coupling h s close to h max at least at the 1% level. The GW signals do not reach 
the sensitivity curve for eLISA, even for the strongest phase transitions. 

We find that the signal from bubble collisions is much weaker than that from turbu- 
lence. As an example, we plot the case g s = 2, fi s = in Fig. 7. The spectrum from 
bubble collisions peaks at a higher frequency. Therefore, bubble collisions will cause a 
secondary peak in the total GW spectrum. However, this peak cannot be observed by 
LISA. For comparison, we show in Fig. 6 the peak signals from turbulence and bubble 
collisions for some values of h s which give a sizeable signal. This result agrees with Ref. 
[22], where only the signal from bubble collisions was considered. 

Bubble collisions may produce a larger signal if the bubble wall can run away. In 
the present model, the effective potential has a barrier at zero temperature for large h s . 
Indeed, the maximum = becomes a false minimum for strong couplings. In this 
context, it is important to ask whether the bubble walls can run away. We have used 
Eqs. (66)-(68) to check for the possibility of runaway walls 10 . We show the result in Fig. 
8 for the case g s = 2. In order to obtain runaway walls, it is necessary to get very close 
to /i max (within a fraction Ah/h ~ 10~ 4 ), and therefore requires a significant fine tuning 
(furthermore, for h s so close to h max the duration of the phase transition quickly becomes 
At 3> H~ v ). Notice that, in Figs. 5-7, the values of h s are not so finely tuned (cf. the 
values of h s for the dots in Fig. 6 and the values reached in Fig. 8). 

5.2 The MSSM 

An interesting example of adding bosons in order to increase the strength of the phase 
transition is the Minimal Supersymmetric Standard Model (MSSM). This model has been 
considered for several years, either in the subject of electroweak baryogenesis (see, e.g., 
[62]) or in that of GW generation [18]. The MSSM contains two complex Higgs doublets 

9 For references on the sensitivity curves see section 6. 

10 We only considered the case \i = 0, which gives more strongly first-order phase transitions. 
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Figure 5: The energy density (top) and frequency (bottom) at the peak of the GW 
spectrum from turbulence, as a function of h s for g s = 2 (rightmost curves) and g s = 12 
(leftmost curves), with /i s = (solid lines), 100 GeV (dashed lines), and 200 GeV (dotted 
lines). Horizontal dotted lines indicate the approximate values corresponding to the peak 
sensitivity of LISA and eLISA, / ~ 1 mHz, h 2 n GW «5x 10~ 12 for LISA, h 2 n GW w 
2 x 10~ 10 for eLISA. 
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Figure 6: The peak GW signal as a function of the peak frequency for the models consid- 
ered in Fig. 5, together with the sensitivity curves for LISA and eLISA. The dots on the 
red lines correspond (from right to left) to h s = 1.94, 1.943, 1.945, and 1.947. The lower 
solid red line is the signal from bubble collisions. 

Hi and H2. We define the vacuum expectation values v\ = (H®) and t>2 = (H®)- ft i s 
customary to simplify the problem by considering the limit in which the CP-odd Higgs 
mass is large (m^ ^> mz)- In this limit the low energy theory contains a single Higgs 
doublet $, and the masses and couplings depend on tan/3 = v 2 /vi. Thus, calling (j)/y2 
the background of the real neutral component of the tree-level potential is of the form 
(56), with the quartic coupling given by A = [g 2 + g' 2 ) cos 2 (2/3)/8. Therefore, the tree- 
level Higgs mass is bounded by m 2 H < m 2 z . However, this tree-level relation is spoiled by 
radiative corrections (see e.g. [63]) and we shall consider itlh as a free parameter. The 
relevant SM field-dependent masses are those of the gauge bosons, m^(0) = g 2 2 /4 = 
h 2 w 2 , m%{4>) = {g 2 + # /2 )0 2 / 4 = /i|0 2 , and top quark, m 2 (0) = h 2 sin 2 /30 2 /2 = h 2 (p 2 , 
where h t is the Yukawa coupling to H®. We shall work in the limit in which the left 
handed stop is heavy {mq > 500 GeV). In this case, the one-loop correction to the SM 
is dominated by the right-handed top squark contribution, with a field-dependent mass 
given by m 2 (<p) ~ mfj + h 2 <p 2 , where 

^ = 0.15/i| cos 2/3 + h 2 (l - A 2 /m 2 Q ^ , (69) 

mfj and mg are soft breaking parameters, and A t is the stop mixing parameter. If the 
mass of the right-handed stop is of the order of the top mass or below, the one-loop 
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Figure 7: The peak energy density from turbulence (solid) and bubble collisions (dashed), 
as functions of h s for g s = 2 and /j, s = 0. 
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Figure 8: The pressure difference p_—p + (dashed) and the mean field value p_—p + (solid), 
as functions of h s for g s = 2 and /i s = 0. The curves are plotted up to h s = 1.9503. The 
total force per unit area p_ — p + becomes positive at h s 1.9497. 
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effective potential admits the high-temperature expansion [62] 



V T (0) = D(T 2 - T 2 ) 2 - T E SM ^ + 6 



12n 



+ 4 ' 



(70) 



where D = m\j (8v 2 ) + 5h 2 w /l2 + 5/i|/24 + h 2 /2 [18], T 2 = m 2 H /(AD), E SM is the 
cubic-term coefficient in the high-temperature expansion for the SM effective potential, 
Esm ~ (2/i^, + h 3 z ) /6n, and M. 2 (0) = m~ (0) +II^(T). The thermal mass is given by 



n t -(T) 
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where t/ s is the strong gauge coupling. We shall set A t = for simplicity in the numerical 
calculation. The parameter T gives the temperature at which the barrier between minima 
of the one-loop effective potential disappears. The phase transition strength is maximized 
for negative values of the soft mass squared mfj ps — IIj (T) [64], for which the contribution 
of the term M."~ in (70) is of the form —E M ssmT4> 3 , with a coefficient E MSS m oc h~ 
that may be one order of magnitude larger than that of the SM. However, such large 
negative values of mfj may induce the presence of color breaking minima at zero or 
finite temperature [65]. In order to avoid the presence of color-breaking minima, we only 
consider values of mfj for which mfj + IT^ (T ) > [18]. 

Nevertheless, the two-loop corrections can make the phase transition strongly first- 
order even for m\j ~ [66]. The most important two-loop corrections are of the form 
4> 2 log <fi and are induced by the SM weak gauge bosons, as well as by stop and gluon loops 
[66, 67]. In the case of a heavy left-handed stop we have [62] 
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where the scale A H depends on the finite corrections and is of order 100 GeV. Following 
[18], we will set A H = 100 GeV for the numerical computation, given the slight logarithmic 
dependence of V2 on A#. 

In the high-temperature approximation, the friction coefficients (62) and (64) become 
[25, 37, 55, 56], 
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(73) 
(74) 



bosons 



where \% = 2 for fermions and \i = m i (0) /T for bosons, a is the surface tension of the 
bubble wall, m 2 D ~ h 2 T 2 is the Debye mass squared, and L w is the width of the bubble 
wall, L w PS 2 /o\ The main contributions to the friction come from the top and the stop. 
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Figure 9: The peak of the GW spectrum as a function of the stop mass, for three values 
of tan (3. 

We consider a range of values of my (corresponding to stop masses in the range 
^stop ~ 130 — 180 GeV) which allow the high-temperature expansion (70) and avoid color- 
breaking minima. We find that the bubbles grow as deflagrations, with wall velocities 
v w ~ 0.4 — 0.5 at the percolation time (slightly higher than at the onset of nucleation 
[37]). The shock-front velocity is v s ^ ~ 0.58. In Ref. [68] a high friction (5 to 10 times 
larger than in an SM-like situation) was obtained using a linear extrapolation of the 
friction from a previous calculation. This gives wall velocities one order of magnitude 
smaller than ours, for which the intensity of GWs would be smaller. 

Figure 9 shows the peak frequency and intensity of the GWs as a function of the stop 
mass for some values of tan/3. The results are quite insensitive to the value of tan/3 for 
tan/3 ~ 1 or higher. For smaller values of tan/3, the phase transition is weaker and the 
intensity of GWs decreases. The results do not change significantly with m stop either. 
The characteristic frequency is f p ~ 20 — 40 mHz. The intensity of the waves is quite low, 
h 2 Q p < 10~ 18 , several orders of magnitude below LISA's sensitivity. This is essentially 
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due to the fact that the coupling of the stop to the Higgs (hf 0.7) is relatively low (cf. 
Fig. 5). The use of a negative mass squared mfj and the two loop correction do not make 
the phase transition strong enough to produce a significant GW signal. 

The results should improve in the Next to Minimal Supersymmetric Standard Model 
(NMSSM), which consists of adding a gauge singlet to the MSSM [18, 69]. A singlet 
extension of the MSSM (the nearly Minimal Supersymmetric Standard Model, nMSSM) 
was considered in Ref. [20], finding that the GW signal is always too low to be observed 
by LISA or BBO. Our results for singlet extensions are more optimistic. We expect in 
this case a similar result to that of adding a singlet scalar to the SM, which we considered 
in section 5.1. The essential difference with the work [20] seems to be the fact that we 
considered the largest bubbles instead of those corresponding to the maximum of the 
volume distribution. As we have seen, this gives an enhancement factor \og(/3/H) in the 
bubble radius. A larger radius decreases the peak frequency as ~ l/Rb but increases the 
GW intensity quadratically 11 . 

Furthermore, in the NMSSM there may be cubic terms, which arise as supersymmetry- 
breaking soft terms. In such a case, the strength of the phase transition is dominated 
by the cubic terms in the tree-level potential, and it is not necessary to rely on loop 
corrections or to consider a light stop. As already mentioned, tree-level effects may lead 
to runaway walls and a larger signal from bubble collisions. 

5.3 Strongly coupled extra fermions 

Extra fermions strongly coupled to the Higgs field can also make the phase transition 
strongly first-order [70] . Strongly coupled fermions, however, make the vacuum unstable. 
This problem can be solved by adding heavy bosons with the same couplings but with 
a large ^-independent mass term, so that they are decoupled from the dynamics at T ~ 
100 GeV. The model can be considered as a particular realization of split supersymmetry, 
where the standard relations between the Yukawa and gauge couplings are not fulfilled. 
In the simplest case, only gf — 12 d.o.f. are coupled to the SM Higgs, with degenerate 
eigenvalues of the form rrij (0) = fi 2 + h 2 (p 2 . Perturbativity requires hf < 3.5. The 
bosonic stabilizing fields have the same number of d.o.f., and a dispersion relation m 2 (0) = 
fi 2 s + h 2 (f) 2 , with h s = hf. For simplicity, Ii s = is assumed. Following [70], we shall set 
ji s to the maximum value consistent with stability, 



In Fig. 10 we have plotted the peak frequency and peak intensity of GWs as a function 
of hf, for several values of \if. Notice that, for high values of the Yukawa coupling hf, this 
model gives mHz frequencies and a GW signal h 2 VL p ~ 10~ 15 , stronger than the MSSM. 
However, the signal is still below LISA's sensitivity. The problem with the extra fermions 
is that, compared to the case of bosons, larger values of the coupling hf are needed to 

11 Notice also that reaching LISA's sensitivity in section 5.1 required a certain fine tuning of the pa- 
rameters. The analysis of Ref. [20] uses randomly chosen values of the parameters, which may not enter 
the fine tuning region. 
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obtain a strongly first-order phase transition. Larger values of hf cause a larger friction 
coefficient. As a consequence, the wall velocity is smaller than in models with extra bosons 
(for a phase transition of the same strength). We find velocities v w < 0.2, and as small 
as v w = 0.05 for strongly first-order phase transitions. This makes this model interesting 
for baryogenesis, since the generated baryon asymmetry peaks for w ffi C 1 [71], but not 
for GW generation. 
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Figure 10: The peak of the GW spectrum as a function of hf for several values of /if. 

6 Detect ability of electroweak gravitational waves: 
LISA and beyond 

In this section we shall compare the results for the models we have considered for the 
electroweak phase transition, and we shall discuss the detectability of the predicted grav- 



32 



itational waves. We show in Fig. 11 some representative curves from each of the models, 
together with the projected sensitivities of several detectors. For comparison, we also 
show other sources of a stochastic GW background, such as galactic and extragalactic 
binaries [72] and inflation. The CMB and large-scale structure constrain the scale of in- 
flation to be below 3.4 x 10 16 GeV, fixing the largest signal expected from inflation [10] to 
/i 2 f2 GW ~ 10~ 14 . Interestingly, for the models we considered, GW signals which are high 
enough to be detected by LISA, are separated in frequency from the noise of white dwarf 
binaries. 

In general, a GW signal of electroweak scale origin lies far away from the sensitivities 
of ground-based detectors such as LIGO or its successors Advanced LIGO, LIGO III, 
which peak at / ~ 100 Hz. Therefore, we shall consider spaceborne detectors. The 
sensitivity curves in Fig. 11 are approximate. The sensitivity for eLISA (upper blue line) 
was calculated using the analytical approximation from Ref. [3]. The other sensitivities 
were calculated from the specifications of the detectors, following the method described in 
Ref. [73]. Specifications for LISA (lower blue line) can be found in [2, 73], specifications 
for BBO (upper purple curve) can be found in [4] and specifications for DECIGO (upper 
orange curve) can be found in [5, 6]. Being composed of several LISA type detectors, the 
latter two will be able to make a correlation analysis between two independent detectors. 
The correlation analysis is expected to increase the sensitivity of BBO to a stochastic 
background by four orders of magnitude [10, 32, 74] (lower purple curve). On the other 
hand, the ultimate sensitivity of DECIGO is estimated to be /i 2 f2 GW ~ 10~ 20 around 
0.1 Hz [5, 10] (lower orange curve). 

The predicted signals for the different models are shown in black in Fig. 11. For all 
the models we considered, the parameters which give frequencies at the sensitivity peak 
of LISA (/ ~ 1 mHz) give intensities a few orders of magnitude below the peak sensitivity 
/i 2 fi G w ~ 10~ 12 . We see that LISA's sensitivity curve is instead achieved at characteristic 
frequencies f p ~ 10~ 4 Hz, by somewhat extreme models, namely, those with extra scalars 
with quite strong couplings to the Higgs. Subsequent detectors like BBO or the Japanese 
DECIGO will have a sensitivity peak about two orders of magnitude below that of LISA, 
/i 2 fi G w ~ 10~ 14 -10~ 13 . However, this peak sensitivity will be for a frequency / ~ 0.1 Hz — 
1 Hz, far away from electroweak GW signals of that intensity. As can be seen in the figure, 
neither BBO nor DECIGO will, in principle, improve LISA's possibility of detecting 
a GW signal from the electroweak phase transition. Nevertheless, after a correlation 
analysis, BBO will possibly be able to detect electroweak GWs for a wider range of SM 
extensions, e.g., extra bosons with moderate couplings or strongly coupled extra fermions. 
The detection would be further improved by the ultimate sensitivity of DECIGO. The 
latter seems to be the only possibility for detecting electroweak gravitational waves in the 
case of the MSSM. In the case of the NMSSM we expect a signal similar to that of the 
SM with an extra singlet. 

7 Conclusions 

We have calculated the intensity and characteristic frequency of gravitational radiation 
generated in the electroweak phase transition. We have considered several extensions of 
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Figure 11: The predicted value of Oqw at the peak of the spectrum as a function of 
the peak frequency, for different models, together with the noise from other stochastic 
sources, and the sensitivities of several spaceborne laser interferometer gravitational wave 
observatories. The upper blue line is the sensitivity curve for eLISA, the lower blue 
line is the one for LISA, the upper purple line is for BBO, the lower purple line for BBO 
correlated, the upper orange line for DECIGO, and the lower orange curve for the ultimate 
sensitivity of DECIGO. The dotted green line corresponds to the signal from white dwarf 
binaries (WD), and the dashed green line is the maximum signal expected from inflation. 
Dotted black curves correspond to SM extensions with extra scalars. From left to right, 
we have g s = 2 d.o.f. with invariant mass fi s = 200 GeV and fi s = 0, and g s = 12 
with fi s = 100 GeV. The intensity of GWs increases with the coupling h s to the Higgs. 
The solid black curve corresponds to the MSSM for tan = 1. The intensity of GWs is 
higher for lower values of the stop mass. Dashed black lines correspond to extensions with 
strongly coupled fermions with gf = 12 and /i/ = (leftmost curve) and /i/ = 300 GeV 
(rightmost curve). In these curves, the GW signal increases with hf. 
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the Standard Model which provide strongly first-order phase transitions, and we have 
discussed the detectability of these models by planned spaceborne gravitational wave 
detectors. 

We have improved the treatment of previous works on the dynamics of the phase tran- 
sition by including in the calculation the hydrodynamics and microphysics of the bubble 
walls. Most works on GWs assume that the bubble walls propagate as Jouguet detona- 
tions. In contrast, we have determined, as a function of the temperature, whether the 
walls propagate as subsonic or supersonic deflagrations, or as weak detonations. We have 
also taken into account the possibility that, instead of reaching a stationary state, the 
walls run away. In order to determine the hydrodynamic solution we have estimated the 
friction for each model, using approximations derived in Refs. [25, 37]. These approxima- 
tions do not depend on details of the specific model and have the correct dependence on 
the relevant parameters (e.g. the couplings of the extra particles with the Higgs). Thus, 
our wall velocity depends on the friction coefficient as well as on the thermodynamic 
parameters. 

Furthermore, we have numerically evolved the phase transition from the nucleation of 
the first bubbles until the time of bubble percolation, taking into account the variation 
of the nucleation rate and of the wall velocity with temperature. We have also taken into 
account the fact that the nucleation rate is suppressed in the regions that are reheated by 
shock fronts. We accomplished this in a simple way by considering the fraction of volume 
occupied by "shock front bubbles" . 

The evolution of the phase transition was considered in some detail in Refs. [20] and 
[22] . The latter studied an extension of the SM with extra scalars. However, only bubble 
collisions were considered as a source of GWs. As we have seen, the signal from bubble 
collisions is generally much lower than that from turbulence (at least for phase transitions 
which do not have runaway bubble walls). The work of Ref. [20], on the other hand, con- 
sidered both signals from bubble collisions and turbulence, and models which, in principle, 
may give stronger phase transitions. However, the relevant size scale of the turbulence 
was assumed to correspond to the bubbles which maximize the volume distribution at 
the percolation time. In contrast, we have argued that the relevant wavelength is given 
by the size of the largest bubbles. We obtain a higher signal, since the intensity of the 
GWs is higher for larger bubbles (pew ~ Rf)- ^ s we nave seen, the size of the largest 
bubbles can be approximated by ps 3t>„,/3 -1 log(/3/H), and there is an enhancement 
\og((3/H) with respect to the size corresponding to maximum volume. It is clear that 
further investigation is needed in order to determine the spectrum of turbulence in the 
presence of several stirring scales (in the case of a phase transition, a continuum of bubble 
sizes). 

For most of the models and parameters, the gravitational wave signal from the elec- 
troweak phase transition seems to be rather weak to be detected by LISA. Nevertheless, 
extensions with scalar singlets which are strongly coupled to the Higgs give considerably 
strong phase transitions, which produce GWs with intensities as high as h 2 il GW ~ 10~ 8 
for frequencies / < 10~ 4 Hz. These models give a signal detectable by LISA, although 
some fine tuning of the parameters (below the 1% level) is required. Taking into account 
that the sensitivity curves are only approximate, and that current calculations of GW 
generation and phase transition dynamics may have large errors, this fine tuning may 
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be relaxed in the future. The extension of the SM with strongly coupled fermions gives 
weaker signals, which could be detected after a correlation analysis from BBO. For the 
case of the MSSM, the ultimate sensitivity of DECIGO would be needed to detect GWs 
from the electroweak phase transition. 

Interestingly, the model with extra fermions gives a larger signal than the MSSM, even 
though the wall velocity is smaller. This confirms the importance of taking into account 
the complete dynamics of the phase transition. To begin with, the wall velocity is not 
directly related to the strength of the phase transition. For instance, two models may have 
the same amount of supercooling and quite different friction, thus giving different wall 
velocities. Most importantly, the GW intensity further depends on the size of the largest 
bubbles, which is rather unpredictable without a careful analysis, due to the nontrivial 
dynamics of nucleation and reheating. Thus, higher wall velocities do not always guarantee 
larger bubble radii. 
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